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6.1 


ABSTRACT 

A  new  two-dimensional  data-adaptive  algorithm  utilizing  the  iterative  Toeplitz 
approximation  method  is  presented.  This  algorithm  provides  a  robust  and  efficient 
means  for  accurate  estimation  of  2-D  autoregressive  parameters  representing  spatially 
variant  data  arrays.  Its  convergence  performance  is  comparable  to  that  of  the  2-D 
Recursive  Least  Squares  (RLS)  algorithm  but  is  computationally  more  efficient.  The 
savings  in  computation  is  realized  by  reducing  the  size  of  the  matrix  to  be  inverted 
when  solving  the  AR  model  normal  equations.  The  ability  of  the  algorithm  to  accu- 
rately estimate  the  model  parameters  using  very  small  data  sets  and  various  window- 
ing schemes  are  evaluated.  Spectral  estimates  are  compared  for  quarter-plane  (QP), 
nonsymmetric  half-plane  (NSHP)  and  combined-quadrant  (CQ)  regions  of  support. 
Additionally,  the  algorithm  is  tested  in  noise  cancellation  and  line  enhancement  ap- 
plications using  image  arrays.  This  algorithm  may  be  implemented  for  data-adaptive 
image  processing  or  coding  and  for  applications  requiring  2-D  spectral  estimation. 
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I.  INTRODUCTION 

A.      PARAMETER  ESTIMATION  OF  DATA-ADAPTIVE  ARRAYS 

Two-dimensional  (2-D)  autoregressive  (AR)  modeling  provides  a  means  for  the 
spectral  estimation  of  signals  received  by  an  array  of  spatially  distributed  sensors. 
Additionally,  the  2-D  AR  model  parameters  may  be  used  to  represent  the  original 
signal  in  a  compact  manner  resulting  in  compression  of  the  data  to  be  processed  or 
transmitted.  This  is  particularly  useful  in  many  image  processing  applications.  The 
traditional  least  squares  or  Wiener  filtering  based  2-D  AR  parameter  estimation  [Ref. 
1]  provide  accurate  spectral  estimates  but  require  the  the  inversion  of  a  relatively  large 
autocorrelation  matrix.  This  is  undesirable  since  matrix  inversion  is  computationally 
expensive.  Iterative  methods,  using  a  Toeplitz  approximation  algorithm,  have  been 
suggested  [Ref.  2,  3]  as  an  alternative  to  the  direct  Wiener  filter  solution.  The  itera- 
tive Toeplitz  approximation  method  has  proven  to  be  an  effective  means  to  estimate 
the  parameters  of  fixed  data  arrays  without  having  to  invert  the  correlation  matrix. 
This  method  is  particularly  useful  when  data  arrays  of  limited  size  must  be  processed. 

The  results  obtained  using  iterative  methods  to  obtain  the  AR  parameters  of 
fixed  data  arrays  suggest  that  these  methods  may  be  used  to  some  advantage  when 
estimating  parameters  adaptively.  Adaptive  parameter  estimation  is  necessary  in 
many  digital  signal  processing  applications  where  the  data  is  non-stationary.  In  these 
applications,  real-time  or  near  real-time  processing  is  often  desirable.  This  requires 
that  any  algorithm  used  must  be  computationally  efficient  and  be  capable  of  providing 
accurate  estimates  obtained  from  a  minimum  of  data.  This  thesis  addresses  the 
implementation  of  the  iterative  Toeplitz  approximation  method  for  2-D  parameter 


estimation  of  non-stationary  data.  It  is  shown  that  this  algorithm  can  be  efficiently 
implemented  to  obtain  accurate  estimates  from  very  small  data  arrays. 

B.  EXPLANATION  OF  NOTATION 

All  rectangular  matrices  are  denoted  by  a  boldface,  captial  letter,  e.g.,  R.  Col- 
umn vectors  are  designated  by  a  boldface  lowercase  letter,  for  example:  x.  Occasion- 
ally, a  vector  created  temporarily  for  the  development  of  a  mathematical  expression 
will  be  represented  by  a  lowercase,  boldface,  greek  letter.  An  example  of  this  would 
be  the  vector  a .  A  special  operator  matrix  will  be  represented  by  a  calligraphic  cap- 
ital letter,  such  as  T.  Scalar  values  are  generally  represented  non-boldface  lowercase 
arabic  or  greek  letters,  e.g.,  a  or  A.  Captital  letters  such  as  A  or  P  are  used  to  rep- 
resent a  variety  of  values  or  expressions.  These  include  anything  from  a  polynomial 
to  the  dimension  of  a  matrix.  The  symbol  T  is  reserved  to  indicate  the  transpose  of 
a  vector  or  matrix. 

C.  THESIS  OUTLINE 

The  following  describes  the  organization  of  the  remainder  this  thesis.  Chapter 
II  serves  as  an  introduction  to  the  problem  of  2-D  AR  modeling  and  parameter  esti- 
mation. Particular  attention  is  given  to  the  formulation  of  the  normal  equations  for 
quarter-plane  (QP)  and  nonsymmetric  half-plane  (NSHP)  regions  of  support.  This 
chapter  concludes  with  a  description  of  the  least  squares  algorithm  and  direct  in- 
version method  of  solving  the  normal  equations.  This  serves  as  the  foundation  for 
Chapter  III  which  extends  the  concept  of  solving  normal  equations  using  the  itera- 
tive Toeplitz  approximation  method  to  derive  parameter  estimates  referred  to  as  the 
fixed  data  method.  Chapter  IV  proceeds  with  the  development  of  the  data-adaptive 
algorithm  using  Toeplitz  approximation.  Two-dimensional  AR  spectral  estimation 
is  discussed  in  Chapter  V,  along  with  a  comparison  of  results  using  both  fixed  and 


adaptive  iterative  methods  for  spectral  estimation.  Additionally,  a  discussion  of  the 
use  of  the  data-adaptive  algorithm  for  noise  cancellation  and  line  enhancement  of 
image  arrays  is  provided  in  Chapter  VI,  with  experimental  results.  Conclusions  and 
recommendations  for  future  work  are  found  in  Chapter  VII. 


II.  2-D  AUTOREGRESSIVE  MODELING 

A.      OVERVIEW 

This  procedure  assumes  a  stationary  (homogeneous)  random  process  x[ni,n2] 
that  is  the  response  of  an  AR  model  excited  by  a  white  noise  input  u>[ni,n2]  having 
a  variance  a1 .  The  AR  model  as  shown  in  Figure  2.1  may  be  described  as  an  all  pole 


w(n1,n2) 


H{u>i,u>2) 


x(ni,n2) 


Figure  2.1:  AR  model  excited  by  white  noise. 


filter  with  the  transfer  function 


1 


where  A  is  the  region  of  support  over  which  the  parameters  a(ki,k2)  are  non-zero. 
The  difference  equation  for  the  system  that  generates  x[ni,n2]  can  be  written  as 


x[ni,n2]  =  -JUIa(i'i):r[ni  -i,n*-  J]  +  tu[ni,n2]. 
(•J)  €A 


(2.2) 


A  linear  set  of  equations  for  the  filter  coefficients  a(i,j)  may  be  formed  by  multiplying 
both  sides  of  (2.2)  by  x\n\  —  h,n2  —  l2]  and  computing  the  the  statistical  expectation 
of  the  resulting  expression  [Ref.  4].  This  leads  to  the  following  expression  called  the 
normal  equation 


^x[/i,/2]  =  -EEa(^i)^[/i-i,/2-i], 

(M)  eA 


(2.3) 


for  /i,/2  >  0.    The  coefficients  a(i,j)  can  be  derived  from  this  normal  equation.    It 
should  be  noted  that  the  structure  of  the  normal  equation  depends  on  the  region  of 
support  A.   A  region  of  support  may  have  any  shape,  but  the  most  commonly  used 
are  the  quarter-plane  and  nonsymmetric  half-plane. 
1.      Quarter-Plane  Support 

The  most  straightforward  region  of  support  is  the  quarter-plane.  A  region 
A  is  considered  to  have  quarter-plane  (QP)  support  when  a(i,j)  has  non-zero  values 
in  one  quadrant  only  as  shown  in  Figure  2.2.   For  this  case  the  normal  equation  has 


n, 
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Figure  2.2:  The  quarter-plane  region  of  support. 


the  form 


Rx[h,h]=  £  x>(i\j)jy/i-t,7a-j],  (i,i)^ (o,o) 

t=0     j=0 


(2.4) 


where/!  =  1,2,  ...,Pi  —  1  and  l2  =  1,2, . . .  ,P2— 1  with  Pi  and  P2  being  the  dimensions 
of  A.  If  it  is  assumed  that  a(0,0)  =  1,  we  may  express  (2.4)  in  block  matrix  form  as 


Ro 

R-! 

R-2 

Ri 

Ro 

R-i 

R2 

Ri 

Ro 

.Rpj-i     Rf,-2     Ra-: 


R-Pi+i 

R-Pi+2 
R-Pl+3 

Ro 


ao 

r^°n 

»i 

0 

a2 

= 

0 

.aPi-i. 

0 

(2.5) 


Each  block  Ra/  of  this  matrix  is  given  by 


Rm  =  R_»^  = 


-M 


R(M,0) 

R{M,1) 


R(M,-1) 

R(M,0) 


i?(M,-P2  +  l) 
R(M,-P2  +  2) 


(2.6) 


R(M,P7-1)     R(M,P2-2)     •••  R(M,0) 

while  the  coefficients  are  a*f  =  [a(M, 0)  a(M,  1)  a(M, 2) . .  .a(M,  P2  —  1)]T  and  the 
error  variance  vector  is  S^  =  \a\  0 . . .  0]T.  The  blocks  Raj  along  the  diagonals  of  the 
autocorrelation  matrix  R  are  equal  and  the  diagonal  elements  of  R^  are  equal,  thus 
R  is  Toeplitz  block- Toeplitz. 

2.      Nonsymmetric  Half-Plane  Support 

A  region  A  with  non-zero  a(i,j)  as  shown  in  Figure  2.3  is  considered  to 
have  nonsymmetric  half-plane  (NSHP)  support.    The  normal  equation  for  this  case 


n? 

i 

A 

i 

' 

i 

Figure  2.3:  The  nonsymmetric  half-plane  region  of  support. 


must  be  modified  to  account  for  non-zero  a(0,j)  in  the  fourth  quadrant.   This  may 

be  written  as 

L2+F2-1  P1-1P2-1 

R*Vu  h]  =      £     a(0,  j)Rx[lu  l2  -  j]  +  J2    £  «(hJ)R*ih  ~  i,  h  -  j)  ,        (2.7) 
j=l  t=0    >=o 


where  (z,j)^0,  h  =  0, 1,2, . . . ,  Pi  -  1  and 
'  0,1,2,...,Z,2  +  P2-1, 


h  = 


for  h  =  0 
L2 ,  L2  +  1 ,  L2  +  2, . . . ,  L2  +  P2  -  1 ,     for  /j  ^  0 


The  dimensions  of  A  are  given  by  P\  and  P2  and  L2  is  a  negative  number  as  defined 
in  Figure  2.3.  As  in  the  quarter-plane  case,  if  we  let  a(0,0)  =  1  we  may  write  (2.7) 
in  block  matrix  form  as 
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Ri 
R2 


R-i 
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R-2 

R-i 
Ro 


R-Pi+i 

R-Pi+2 
R-Pi+3 


ao 
ai 

a2 

= 

■S(oy 

0 
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.apj-i. 

0 

(2.8) 


.Rpi-i    Rp!-2    Rpi-3    •••       Ro 

The  blocks  Ro,  R^  and  Rm  each  have  different  forms.  The  diagonal  block 
Ro  has  the  form 


Ro  = 


#(0,0) 

#(0,1) 


#(0,-1) 
#(0,0) 


#(0,-L2-P2  +  l) 
#(0,-Z,2-P2  +  2) 


#(0,/2  +  P2-l)     #(0,L2  +  P2-2)     •••  #(0,0) 

the  blocks  Raj  =  R^  in  the  first  row  and  column  may  be  written  as 


#(M,-L2) 
#(M,-L2  +  1) 


#(M,-L2-1) 
#(M,-L2) 


#(M,-L2-P2  +  1) 
R{M,-L2-P2  +  2) 

#(M,0) 


#(M,-L2  +  P2-1)     #(M,  -L2  +  P2  -  2)     ••• 
and  the  remaining  blocks  Rm  are  given  by  (2.6). 

The  model  parameter  vectors  are  given  by 

ao  =  [a(0, 0)  a(0, 1)  a(0, 2) . . .  a(0,  L2  +  P2  -  1)]' 
and 


(2.9) 


,(2.10) 


aM  =  [a(M,  L2)  a(M,  L2  +  l)  a{M,  L2  +  2) . . .  a{M,  L2  +  P2  -  l)}: 


for  M  =  1, 2, . . . ,  Pi  - 1.  The  error  variance  vector  £<°>  =  [<r£2  0 . . .  0]T.  Except  for  the 
upper  diagonal  block  (2.9)  and  top  and  left  most  blocks  (2.10),  the  NSHP  autocorre- 
lation matrix  is  block-toeplitz  with  each  block  being  toeplitz  as  well.  Quarter-Plane 
support  can  be  considered  a  special  case  of  NSHP  support  with  L2  =  0. 

B.     SOLVING  THE  NORMAL  EQUATIONS  BY  DIRECT  INVERSION 

The  normal  equations  also  arise  in  the  2-D  linear  prediction  problem.  When 
solving  the  normal  equations,  the  objective  is  to  find  the  parameters  which  mini- 
mize the  prediction  error  [Ref.  5].  It  is  frequently  more  convenient  to  describe  the 
formulation  of  the  normal  equations  in  terms  of  linear  prediction. 

In  the  2-D  least  squares  problem,  the  error  between  the  random  process  x\ri\,  712] 
and  its  estimate  a:  [711,7*2]  is  given  by 

e[7ii,n2]  =  z[ni,n2]  —  x[ni,n2]  (2.11) 

or  in  expanded  form  as 

e[nx,n2]  =  x[nu n2]  +  £  J2 a(ij)x[n!  -  i,n2  -j],  (ij)  ^  (0,0).  (2.12) 

(•Vj)  6^4 

The  objective  of  the  least  squares  method  is  to  minimize  the  squared  errors  from  a 
particular  set  of  these  terms  [Ref.  4].  If  we  let  a(0,0)  =  1,  we  can  express  (2.12)  as 

e[nun2]  =  ]T  ^2a{iJ)x[ni  -  i,n2  -  j]  (2.13) 

which  is  compactly  represented  in  vector  notation  as  e  =  Xa.  The  error  e[ni,n2] 
is  computed  for  each  position  of  the  filter,  then  normal  equations  are  formed  by 
multiplying  both  sides  of  (2.13)  by  XT  and  applying  the  orthogonality  principle  [Ref. 
5].  This  results  in 

Ra  =  S  (2.14) 
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where  R  =  XTX  and  the  error  variance  vector  S  =  XTe  =  [£{o)T,0, . . .  ,0]T  with 
£(°)  =  [<t£2,  0, . . . ,  0]T.  The  matrix  X  may  be  defined  as 


X  = 


xo     X' 


X0       X! 


XPj-l 


(2.15) 


rp 

and  the  model  parameters  may  be  given  by  a  =  [  1  a'  ]  .  The  parameter  estimates 
are  obtained  from  (2.14)  by  muliplying  both  sides  of  the  equation  by  the  inverse  of 
R  which  may  be  written  in  expanded  form  as 


a'  =  -(XriXr1XAJx0 


(2.16) 


This  is  referred  to  as  the  direct  inversion  method. 

The  preceding  discussion  assumes  that  R  is  Toeplitz  block-Toeplitz.  This  true 
for  the  case  where  the  autocorrelation  method  [Ref.  5]  is  used  to  form  the  correlation 
matrix  or  when  the  theoretical  correlation  is  known.  In  practicality,  however,  param- 
eter estimates  must  be  obtained  from  relatively  small  sets  of  finite  data.  In  this  case, 
it  is  often  desirable  to  compute  R  using  the  covariance  method.  This  results  in  a 
block  autocorrelation  matrix  R  with  non-Toeplitz  blocks  R^  [Ref.  2].  The  rows  of 
the  matrix  X  are  the  elements  of  the  random  process  covered  by  the  filter  mask  for 
each  point  being  estimated  as  the  mask  is  moved  over  the  data.  In  the  covariance 
method  the  filter  mask  is  not  moved  past  the  boundaries  of  the  region  of  support. 
This  means  that  when  QP  support  is  used,  X  will  be  of  dimension  P\Pi  x  P\Pi  where 
Pi  and  P2  are  the  dimensions  of  the  data  array  accessible  to  the  filter  mask.  For 
example,  using  a  3  x  3  (nine  element)  filter  mask  to  form  the  normal  equations  of  an 
8x8  data  array  would  result  in  an  X  matrix  with  dimension  36  x  9.  It  should  be 
noted  that  Px  =  P2  =  6,  since  two  rows  and  two  columns  of  data  cannot  be  estimated 
by  the  filter. 

The  preceding  discussion  provides  the  basis  for  the  subsequent  work  in  this 


thesis.  The  next  chapter  will  addresses  this  same  problem  using  the  iterative  Toeplitz 
approximation  method. 
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III.  ITERATIVE  TOEPLITZ 
APPROXIMATION  ALGORITHM 

The  major  drawback  of  the  least  squares  method  described  in  the  previous 
chapter  comes  from  the  requirement  to  invert  a  relatively  large  autocorrelation  ma- 
trix. This  is  undesireable  since  matrix  inversion  is  computationally  expensive.  In 
this  chapter  we  present  an  iterative  method  for  solving  the  normal  equations  which 
does  not  involve  the  direct  inversion  of  the  correlation  matrix  to  obtain  the  model 
parameters.  This  method  is  then  extended  to  the  case  where  the  covariance  method 
is  used  to  form  the  correlation  matrix  and  the  normal  equations  are  solved  iteratively 
using  Toeplitz  approximation. 

A.     THE  ITERATIVE  SOLUTION  OF  THE  NORMAL  EQUATIONS 

An  alternative  to  direct  inversion  is  to  take  advantage  of  the  near  Toeplitz- 
block- Toeplitz  structure  of  the  autocorrelation  matrix  and  successively  partition  the 
normal  equation  [Ref.  2].  This  partitioning  permits  the  estimation  of  parameters 
while  requiring  the  inversion  of  a  matrix  that  is  significantly  smaller  than  the  original 
autocorrelation  matrix.  The  following  discussion  develops  the  iterative  solution  for 
QP  and  NSHP  regions  of  support. 

1.      Quarter-Plane  Support 

In  order  to  iteratively  solve  the  QP  normal  equations  we  must  begin  by 
dividing  both  sides  of  (2.5)  by  a2.  This  results  in  a  modified  normal  equation  given 
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by 


Ro 

R-i 

R-2 

Ri 

Ro 

R-i 

R2 

Ri 

Ro 

R-A+i 

R-Pi+2 
R-Fi  +3 


'    ag    1 
a'/ 

. 

-£(0)- 

0 
0 

L  aPj  _i  J 

0 

(3-1) 


Rpj_i    Rpj-2    Rpi-3    •••       Ro 
where  £(0)  =  [1  0  0, . . . ,  0]T  and  a?  =  £a<  for  i  =  1,2, . . .  ,i\  -  1.  We  now  partition 

the  normal  equation  (3.1)  as  follows 

(3.2) 


where 


\G1     ha] 

X   rJ 

«i 

a" 

L  aP_j  j 

= 

"*r  i" 
.  0  . 

G1  = 


Ro  R-i  R-2 

Ri         Ro  R-i 

R-2         Ri  Ro 

RPi  -2  RPi  -3  RPi  -4 


R-Pi+2 
R-Pj+3 
R-Pi  +4 


Ro 


(3.3) 


hj ■=  [Rpi-i    Rpj-2    Rpi-3    ■••    Ri]  , 


(3-4) 


and 


—  r  -."^     *"T        ~iit 


o  i  =  I  aj,      aj 


lPi-2l       > 


T\T 


(3.5) 


7  ^[f^7     0T     ...     0TY   .  (3.6) 

We  may  now  form  a  set  of  coupled  iterative  equations  by  defining  an  op- 
erator T\  which  is  given  as 


?i  = 


Gr1    o 

0       Ro-1 


and  using  this  operator  to  premultiply  both  sides  of  (3.2)  to  yield 


°i  =  Gi  hi  _hi  a?k-i]  > 


and 


ifi-i  = 


-luT. 


-Ro-'hja! 
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(3.7) 


(3-8) 


(3.9) 


Equations  (3.8)  and  (3.9)  suggest  an  iterative  solution  and  can  be  written  as 


a 


(*)  _  r--i 


"(*-i)i 


Grxb  i   -hiaSV], 


(3.10) 


and 


a^  =  -Ho-hM*-" 


(3.11) 


which  provide  a  means  to  solve  the  normal  equation  (3.1)  iteratively. 

It  can  be  noted  that  the  submatrix  Gi  is  nearly  the  same  dimension  as 
the  original  correlation  matrix.  This  means  that  the  inversion  of  Gi  provides  little 
advantage,  computationally,  over  the  direct  inversion  method.  It  is  clear,  however, 
that  further  partitioning  of  Gi  will  decrease  the  required  complexity.  Partitioning  of 
Gi  yields 

(3.12) 


where 


a2 


G,  = 


h2      RqJ  L&p_2. 


Ro  R-i  R-  2 

Ri         Ro  R-i 

R2         Ri  Ro 

.Rpj-3  Rpj-4  Rpi-j 


f  2 
0 


R-P1+3 

R-Pi+4 

R-Pi+5 


Ro 


(3.13) 


h2  =  [^1-2    R-Pi-3    Rpi-^ 


Ri]  , 


(3.14) 


and 


—  f  *>ffT     Q'/T 


a2  =  [aj 


J  IT 


JiT      \T 


T      U 
i— 3  J        5 


7   2  =  [£(°)T      0T 


TlT 


oT] 


(3.15) 


(3.16) 


As  before,  both  sides  of  (3.12)  are  premultiplied  by  an  operator  defined  by 


T*  = 


G21 


0 
0       Ro"1 


(3.17) 
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which  results  in 


and 


a  2  =  G1 1  [7  2   -  h2  aPi  _2]  , 


aD.  _0  =  — xtn  n0  Or 


(3.18) 


lPi_2  -  -xt0   n2a2  .  (3.19) 

Equations  (3.18)  and  (3.19)  may  be  solved  iteratively  in  the  same  manner  as  (3.8) 
and  (3.9).  This  solution  is  expressed  as 


and 


«W  =  G2-'h1   -b.  •#?>), 


aPi-2  —  —  "0   "2  a  2 


(3.20) 


(3.21) 


This  partitioning  process  may  be  continued  until  Gp,_i  =  Ro-   This  will 
result  in  a  partitioned  normal  equation  given  by 


[Gp,.!     hp,.!]  fap^i]  _  [7  Pi_! 
[h^.!        Ro    J  LaP_2  J"  L      0 


(3.22) 


where  hp^  =  Ri,  otpj-i  =  a£  and  7  Pi-1  =  S^°\ 

Premultiplication  of  both  sides  of  (3.22)  by  an  operator  Tp^-\  will  result 
in  the  iterative  solution 


„»(*)  =  vc.xe ,0,  _  [R_ia»(*-D  + . . .  +  R_Pi+ia'^)]  , 


and 


a'/""  =  -Ro'IR,^'-'  +  R-,af ->  +  . . .  +  R-^a^"] 
For  this  case  the  operator  matrix  is  given  by 


(3.23) 


(3.24) 


Fft-1  = 


[Ro"1       0 
I    0       R^ 


-1 


(3.25) 
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At  this  point  in  our  development,  we  can  combine  the  preceding  iterative 
solutions  from  the  successive  levels  of  partitioning  into  a  set  of  compact  iterative 
summations 

4'w  =  aT>  -  Ro"1  J:  R-.af*-"  ,  (3.26) 


af  >  =  -R^1  £  Ri-.^"*"11  .  d  *  i)  ,  (3-27) 

i=0 

where  4'(0)  =  R^*0),  af0)  =  -R^R^,  and  j  =  1, 2, . . . ,  P1  -  1.  The  index  of 
iteration  is  k  and  a"  are  the  P<i  x  1  parameter  vectors.  The  solution  of  (3.26)-(3.27) 
requires  0(P12P22^)  multiplications,  where  A"  represents  the  number  of  iterations 
to  converge  to  the  true  parameter  values.  An  additional  P\P%  multiplications  are 
required  to  compute  the  initial  parameter  estimates  a^  for  j  =  1,2, . . .  ,  P2  —  1. 
This  results  in  total  multplications  of  Q{P\PlK  +  PiP22)-  Depending  on  the  number 
of  iterations  required  to  converge  to  the  true  parameters,  this  total  compares  with  the 
algorithms  of  Wax  and  Kalaith  [Ref.  6]  and  Akaike  [Ref.  7]  which  require  (9(P13P22) 
operations.  When  large  arrays  are  used  this  represents  a  noticeable  savings  over  direct 
inversion  which  requires  0(P13P23)  multiplications. 
2.      Nonsymmetric  Half-Plane  Support 

The  derivation  of  the  iterative  NSHP  solution  follows  a  procedure  of  suc- 
cesive  partitioning  similar  to  that  of  QP  support.  However  some  differences  exist  due 
to  the  asymmetry  of  the  autocorrelation  matrix  consisting  of  three  different  types 
of  block  matrices.  This  asymmetry  results  in  a  somewhat  modified  set  of  recursive 
summations  given  by 

a?*>  =  aj(0)  -  Ro"1  ^  iL.a"**-1)  (3.28) 

t=i 

af  >  =  -Ro-Rj'af -"  -  R;1  £  R^a""-'*,  (i  *  j),  (3.29) 

t  =  l 
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where  a£(0)  =  R^#°>,  a;'(0)  =  -R^R.a?0*,  and  j  =  1,2, . . .  ,P,  -  1.  The  index 
of  iteration  is  k  and  a''  are  the  P2  X  1  parameter  vectors.  As  with  QP  support,  The 
solution  of  (3.28)-(3.29)  requires  O^fPfK)  multiplications,  where  K  represents  the 
number  of  iterations  to  converge  to  the  true  parameter  values. 

B.     THE  ITERATIVE  SOLUTION  USING  TOEPLITZ  APPROXIMA- 
TION 

The  discussion  of  the  previous  section  assumes  an  autocorrelation  matrix  with 
Toeplitz-block-Toeplitz  structure.  In  most  cases  the  autocorrelation  matrix  must  be 
formed  from  a  limited,  and  possibly  very  small,  amount  of  sampled  data.  In  this 
situation,  it  is  best  to  use  the  covariance  method  to  estimate  the  correlation  matrix 
to  reduce  the  bias  in  the  estimate.  However,  since  this  method  does  not  have  Toeplitz- 
block-Toeplitz  structure,  it  is  necessary  to  modify  the  iterative  algorithm  described 
above. 

1.      Forward  Iteration 

For  first  quadrant  QP  support,  the  covariance  estimate  R  has  the  form 


Ro,o 

Ro.i 

R<),2 

Ri,o 

Ri.i 

R-1,2 

R-2,0 

R24 

R-2,2 

Ro,Pj-i 
RitPi-i 

R-2,Pi-l 


(3.30) 


Rpi-1,0    Rpi-1,1    Rpi-1,2    •••    Rpi-i,Pi-i 

A  proven  procedure  [Ref.  8]  that  can  be  used  to  form  the  Toeplitz  approx- 
imation is  to  first  average  the  diagonal  blocks  Rj,j,  i  =  j  and  then  form  a  Toeplitz 
T  matrix  from  this  averaged  block  Ravfl  by  averaging  its  diagonal  elements.  The 
diagonal  elements  of  T  can  be  obtained  as  follows 


t(i)  = 


1 


Pj-t-i 


H    RaVp(i+i,j) 


(3.31) 


3=0 


where  i  =  0, 1, . . . ,  P2  —  1. 
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We  may  now  proceed  to  succesively  partition  the  normal  equations  in  the 
manner  outlined  in  the  previous  section.  The  first  partitioning  leads  to  the  iterative 
equations 

a^)  =  Gr1[71   -  hx  a^"1*] ,  (3.32) 

and 

A%  =  -Rp,1-i.P,-ih[a<t-1»  .  (3.33) 

The  diagonal  block  Rpj_i,Pi_i  may  be  written  as 

Rp^i.p,-!  =  T  +  A^-lp,.!  (3.34) 

where  Apj-i.Pj-i  is  the  difference  between  Rpj-i^-i  and  the  Toeplitz  approximation 
T. 

Substituting  (3.34)  into  (3.33)  modifies  the  recursion  to  give 

-f'-Gr'bi   -h.a^.-1'],  (3.35) 

and 

'r(k)     _       rp-lA  «"(*-!)        T-1UT      (k-1)  (O  oc\ 

apj_i  =  -J-     &p1-\,p1-ia.Pl_1    -  1     r^oti         .  {6.6b) 

The  final  expressions  resulting  from  the  successive  partitioning  are  then 

aglj  =  Gp^b  Pl-t   -  hp,.!  a'^"1*]  ,  (3.37) 

and 

a;,(fc)  = -R^hp^xa^V  •  (3.38) 

where  apj_i  =  ao,Gpj_i  =  Ro,o,   and  hpa_i  =  Roj.    The  diagonals  Hjj  may  be 
rewritten  as 

R,,,  =  T  +  Ajj  ,  (3.39) 
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Now,  substituting  (3.39)  into  (3.38)  results  in 


a 


(3.40) 


and 

a/     =  -T     Ai.iaf      ' -T     n1a^_1'.  («*.41) 

Finally,  combining  the  succesive  iterative  equations  results  in  our  solution  to  the 
normal  equations  using  the  forward  iterative  algorithm  for  the  QP  case,  which  takes 
the  form 


a?*>  =  aT  -  T-1A0,o4,(fc-1)  -  T"1  £  RoX 


rt(k-l) 


Pi-l 


//(fc-1) 


t=l 


a 


"(*)  _ 


Pi-i 


=  -T^A^"1'  -  T"1  £  R^a* 


/#(*-!) 


(3-42) 
(3.43) 


i  =  0 


where  a?0*  =  T"1^0),  a"(0)  =  T^R^a?0*  and  j  =  1,2,..., Pi  -  1.   The  forward 
iteration  will  provide  parameters  for  the  first  quadrant  spectral  estimate. 
2.     Backward  Iteration 

A  filter  mask  may  have  quarter-plane  support  in  any  two  quadrants.  A  sec- 
ond quadrant  estimate  exists  for  the  quadrant  adjacent  to  either  side  of  the  quadrant 
designated  to  be  the  first  quadrant.  The  Hermitian  symmetry  and  Toeplitz-block- 
Toeplitz  nature  of  the  autocorrelation  matrix  causes  the  estimates  produced  from 
diagonally  adjacent  quadrants  to  be  identical  [Ref.  9]. 

The  second  quadrant  AR  parameters  b(i,j)  may  be  estimated  from  the  nor- 
mal equations  using  the  backward  iterative  Toeplitz  approximation.  When  6(0, 0)  =  1 
the  second  quadrant  normal  equation  is  given  by 


Ro,o 
Ri.o 

R2,0 


Ro,i 

Ri.i 

R2.1 


Ro,2 
Rl,2 
R2,2 


Rpi-i.o    Rpj-i,i    Rpi-1,2 


Ro.Pi-i 
Ri,Pi-i 

R2,Pi-l 

Rpj-i.Pj-i. 
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'bpj-i" 

0 

bpj-2 

0 

bpj_3 

= 

0 

L    b0    . 

£(0) 

(3.44) 


where  R^  =  R.T.,  £(°)  =  [1  0  0, . . . ,  0]t  and 

bTk  =  (l/*2)[b(k,0),b(k,l),b(k,2),...  ,6(^-1)]. 

It  may  be  noted  that  the  second  quadrant  autocorrelation  matrix  is  identical  to  that 
of  the  first  quadrant.  However,  the  parameter  estimates  b(i,j)  are  not  generally  the 
same.  The  estimation  of  the  backward  parameters  is  computed  iteratively  in  the  same 
manner  as  the  forward  parameters,  but,  the  backward  method  differs  in  the  way  that 
the  normal  equation  is  partitioned.  In  this  case  the  partitioning  begins  at  the  upper 
left  diagonal  block  and  continues  until  the  G  matrix  consists  of  only  the  lower  right 
diagonal  block.  The  Toeplitz  approximation  T  is  identical  to  that  of  first  quadrant 
support. 

The  combined  backward  iterative  equations,  which  are  similar  to  those  for 
the  forward  iteration,  can  be  written  as 

bj*>  =  b<°>  -  T-1Ar1_1,Pl_1b<t-,)  -  T-  £  B*_1A_1_»bj£Jk  (3.45) 

i=l 

bf  =  -T-^.^p^bf-1*  -  T"1  £  Rft-i-iA-^bglV-*  (3-46) 

•=0 

where  b<0)  =  T"1^0),  bf  =  T^R^-ub^  and  j  =  1,2, . . .  ,PX  -  1.  The  first  and 
second  quadrant  parameters  are  used  to  find  the  combined-quadrant  spectral  estimate, 
which  is  discussed  in  Chapter  V. 

3.      Iterative  Method  for  NSHP  Support 

The  iterative  Toeplitz  approximation  method  may  be  used  to  estimate  the 
parameters  of  arrays  with  nonsymmetric  half-plane  support.  As  shown  previously  the 
NSHP  autocorrelation  matrix  has  an  asymmetric  structure  which  can  be  expressed 
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as 


Ro.o 

Ro,i 

R(),2 

Ri,o 

R-1,1 

Rl,2 

R-2,0 

R-2,1 

R-2,2 

Ro,Pi-i 
R-i.Pi-i 


(3.47) 


R-Pi-i.o    Rpi-i.i    R-Pi-1,2    •••    R-Pi-i.Pi-i 

The  dimensions  of  the  top  most  and  left  most  blocks  are  reduced  by  L2  rows  and 
columns,  respectively.  As  in  the  previous  cases  the  iterative  process  begins  with 
forming  the  Toeplitz  approximation  of  Ra«5.  In  the  NSHP  case,  however,  the  upper 
left  main  diagonal  block  is  not  included  in  the  average  becuase  it  does  not  have 
the  same  dimension  as  the  other  main  diagonal  blocks.  Thus,  it  is  necessary  to  use 
the  elements  of  the  upper  left  main  diagonal  block  to  create  a  separate  Toeplitz 
approximation  designated  as  T. 

The  normal  equation  is  successively  partitioned  as  described  in  the  previous 
sections.  At  each  stage  of  the  partitioning,  the  blocks  along  the  main  diagonal  are 
expressed  as  the  sum  of  the  Toeplitz  approximation  and  a  difference  matrix.  The 
iterative  summations  are  given  by 


S"(*)  _.  5»(o) 


P1-1 


T-1  AccST"1'  -  T-1  £  R^aJ 


//(*-!) 


i=i 


"(*)  _ 


=  -T-^a^"1'  -  T^R^a?*"1*  -  T"1  £  R^a," 


s"(*-i) 


P1-1 


//(fc-i) 


(3.48) 
(3.49) 


t=0 


where  a£(0)  =  f-iflW,  «*«  =  T"1^^  and  j  =  1,2, . . . , Px  -  1. 


C.     SUMMARY 

It  is  clear  that  iterative  methods  have  a  computational  advantage  over  the  di- 
rect inversion  method.  This  has  been  shown  for  the  Toeplitz-block- Toeplitz  case  and 
extended  to  the  case  where  Toeplitz  approximation  is  used  to  compensate  for  the 
non-Toeplitz  nature  of  autocorrelation  matrices  formed  using  the  covariance  method. 
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Generally,  the  matrix  that  must  be  inverted  using  the  iterative  method  has  the  di- 
mension of  a  single  block  in  the  block  matrix.  This  corresponds  to  the  filter  order  or 
mask  size.  The  mask  size  will  vary  with  the  order  necessary  for  accurate  AR  mod- 
eling of  the  sampled  data,  but  filter  masks  that  range  from  3  x  3  to  6  x  6  will  be 
sufficient  for  most  applications.  Although  other  methods  may  exhibit  faster  conver- 
gence to  the  true  parameters,  they  do  so  with  greater  computational  complexity.  The 
iterative  methods  provide  an  approximation  much  sooner.  Additionally,  the  iterative 
methods  are  guaranteed  to  converge  when  the  block  matrix  is  symmetric  and  positive 
definite  [Ref.  2].  This  property  is  important  when  extending  the  iterative  Toeplitz 
approximation  algorithm  to  the  data-adaptive  case  which  is  introduced  in  the  next 
chapter. 
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IV.  DATA- ADAPTIVE  ITERATIVE  TOEPLITZ 
APPROXIMATION  ALGORITHM 

A.     OVERVIEW 

The  necessity  for  adaptive  filtering  techniques  arises  when  it  is  desired  to  process 
signals  that  result  from  an  environment  of  unknown  statistics  [Ref.  10].  There  exists  a 
broad  class  of  problems  that  fall  into  this  category,  these  include  such  diverse  fields  as 
sonar,  radar,  image  processing,  seismology  and  communications.  In  general,  adaptive 
filters  provide  a  significant  improvement  in  performance  over  fixed  coefficient  filters 
designed  to  operate  on  data  with  unknown  statistics.  Additionally,  the  use  of  adaptive 
filters  makes  possible  new  signal  processing  capabilties  that  would  not  be  available 
otherwise. 

One  distinct  advantage  associated  with  adaptive  filters  pertains  to  their  ability 
process  data  on-line.  This  is  desirable  for  many  applications  such  as  autoregressive 
spectrum  analysis,  detection  of  a  signal  in  noise,  adaptive  noise  cancellation  and  line 
enhancement. 

Two-dimensional  adaptive  filters  are  used  to  process  data  obtained  from  spa- 
tially variant  data  arrays.  The  most  successful  algorithms  currently  implemented  for 
this  purpose  include  the  2-D  least  mean  square  (LMS)  and  recursive  least  squares 
(RLS)  algorithms.  Both  of  these  algorithms  have  been  derived  from  the  2-D  Wiener 
filter  and  method  of  least  squares  [Ref.  11].  With  this  in  mind,  it  follows  that  the 
iterative  method  for  fixed  data  may  be  successfully  extended  to  operate  on  spatially 
variant  arrays.  The  remainder  of  this  chapter  describes  the  development  of  the  data- 
adaptive  iterative  Toeplitz  approximation  algorithms  and  some  considerations  that 
apply  when  using  this  algorithm. 
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B.     THE  DATA-ADAPTIVE  ALGORITHM 

The  data-adaptive  Toeplitz  approximation  algorithm  is  based  on  the  same  suc- 
cesive  partitioning  described  for  the  fixed  data  case.  The  algorithm  becomes  adaptive 
in  nature  when  the  autocorrelation  matrix  Rn  is  updated  for  each  shift  n  of  the  filter 
mask.  This  yields  a  continuously  updated  set  of  AR  parameter  estimates.  Compu- 
tation of  AR  parameter  estimates  begins  with  the  first  iteration  [Ref.  12].  As  in  the 
fixed  data  case,  only  the  approximated  Toeplitz  matrix  T  is  inverted.  Computation 
time  is  shortened  further  by  the  fact  that  R^  is  computed  only  from  data  covered 
by  the  filter  mask,  which  amounts  to  taking  the  outer  product  of  a  P\Pi  x  1  vector 
where  P1P2  is  equal  to  the  number  of  filter  coefficients.  The  data-adaptive  form  of 
the  combined  iterative  equations  can  be  written  as 

3j(»)  =  aj(o)  _  T-i  A^aJ-D  _  T;1  Y:  Rn.cMaf1"1*  (4.1) 

t=i 

a»W  .  _t;.Aivj  ja»(«-.)  _  T-i  £  Rnj„a;'("-11  (4.2) 

t=0 

where  a'0/(0)  =  T"1^0),  ao"(0)  =  To"1Roil0aS(0)  and  j  =  1, 2, . . . ,  Px  - 1.  These  can  be 
compared  to  (3.42)  and  (3.43)  for  the  fixed  data  case.  The  initial  parameter  estimates 
ao"  are  determined  for  the  initial  point  being  estimated  and  all  subsequent  parameters 
a'-  are  a  function  of  the  parameters  computed  at  the  previous  mask  positions.  The 
backward  equations,  used  for  the  second  quadrant  estimates,  can  be  written  as 

b<n>  =  bf  -  T-1  A^^hj"-1)  -  T;1  J2  Rn,pJ_1,p1_1_,b^:11)_t  (4.3) 

t=l 

b(n)  =  _T-lAnPi_1_.pi_1_.b(n-l)  _  T-l    f:   Rntp1_1_J>p1_1_lb^:V_t  (4.4) 

1=0 

where  b^  =  T"1^,  bo;o)  =  To^Rop,-!,;^  and  j  =  1,2, . . .  ,Pa  -  1. 

It  should  be  noted  that  the  only  difference  between  the  data-adaptive  equations 
and  the  fixed  data  equations  is  the  the  index  of  iteration  which  is  n  instead  of  k 
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and  the  manner  in  which  the  autocorrelation  matrix  is  formulated.  In  the  data 
adaptive  case,  the  iteration  proceeds  with  each  shift  of  the  filter  mask  or  update  of 
the  autocorrelation  matrix  which  has  the  recursive  form 

R*  =  Rn-!  +  XnXJ  .  (4.5) 

The  P1P2  x  1  vector  xn  is  obtained  by  arranging  the  2-D  data  covered  by  the  filter 
mask  at  each  shift  n  =  (ni,n2)  in  a  vector  form.  The  diagonal  blocks  of  Rn  are  used 
to  form  Tn  in  the  same  manner  as  the  fixed  data  case. 

As  with  the  RLS  algorithm  and  other  recursive  implementations  of  the  method 
of  least  squares,  it  is  necessary  to  introduce  a  weighting  or  forgetting  factor  when 
computing  R„  [Ref.  10].  The  forgetting  factor  is  a  scalar  value  that  may  be  designated 
as  /?(n)  where  n  is  the  iteration  number  of  the  point  being  estimated.  The  weighting 
factor  P(n)  has  the  property  that 

0<  0(n)<  1,    n  =  l,2,...,p,  (4.6) 

where  p  is  the  total  number  of  points  being  estimated.  The  purpose  of  fi{n)  is  to 
ensure  that  data  in  the  distant  past  is  'forgotten'  which  will  make  it  possible  for 
the  filter  to  adapt  to  statistical  variations  of  the  observed  data  when  operating  in  a 
non- homogeneous  environment.  A  commonly  used  form  of  the  weighting  factor  is  the 
exponential  weighting  factor  which  is  defined  as 

P(n)  =  \p-n,    n  =  l,2,...,p,  (4.7) 

where  A  is  positive  constant  with  a  value  less  than  1.  The  quantity  1/(1  —  A)  can  be 
considered  as  an  approximate  measure  of  the  memory  of  the  algorithm.  When  A  =  1 
we  have  the  ordinary  method  of  least  squares  which  corresponds  to  infinite  memory. 
By  this  reasoning,  we  may  introduce  a  method  of  exponentially  weighted  least  squares, 

24 


where  fi(n)  is  used  as  a  weighting  factor  in  the  expression  of  the  performance  measure 
[Ref.  10] 

e(n)  =  EA"-"|e(n)|2,  (4.8) 

n=l 

where  e(n)  is  the  error  defined  by  (2.12)  at  mask  position  (ni,n2). 

To  see  how  the  forgetting  factor  is  introduced  in  the  recursion  we  must  refer  to 
the  original  definition  of  the  normal  equation 

Ra  =  €  ,  (4.9) 

where  a  is  the  vector  of  optimum  parameters  which  results  in  the  minimum  value  of 
the  performance  measure  £(n).  The  P1P2  X  P\P<i  autocorrelation  matrix  in  this  case 
is  defined  by 

Rn  =  J2  ^"n*nXj  •  (4.10) 

n=l 

We  may  now  isolate  the  term  of  (4.10)  that  corresponds  to  n  =  p  from  the  rest  of  the 
summation  on  the  right  side  of  the  expression  to  obtain 


Rn  =  A 


n=l 


JT 


+  x„xj.  (4.11) 


The  term  inside  the  brackets  on  the  right  side  of  (4.11)  is  actually  the  sum  of  the  pre- 
viously computed  autocorrelation  matrices  R(n  —  1),  therefore  we  have  the  recursive 
relationship 

Rn  =  AR^  +  xnxj  .  (4.12) 

The  affect  of  0(n)  is  apparent  when  (4.12)  is  expanded  to  give 

Ri     =    Xxxf 

R2    =    Axixf  +  x2x^ 

R3    =    A2xixf  +  Ax2xJ  +  X3X3 


R„    =    A^xjxf  +  An-'x2x^  +  . . .  +  xnxi  (4.13) 
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which  clearly  shows  the  decreasing  contribution  made  by  past  data  vectors  when 
A<  1. 

The  choice  of  a  value  for  A  is  dependent  on  the  statistical  nature  of  the  data 
being  processed.  Generally,  if  the  data  is  random  or  uncorrected  it  is  undesirable  to 
use  the  past  data  to  compute  the  current  estimate,  therefore  a  small  value  for  A  would 
be  the  appropriate  choice.  Conversely,  a  value  of  A  close  to  1  would  be  desirable  for 
highly  correlated  data. 

Since  the  statistical  nature  of  the  data  is  often  not  known  a  priori  when  using 
data-adaptive  filters,  it  is  difficult  to  choose  an  optimum  value  for  the  forgetting 
factor.  A  further  modification  of  the  recursive  equation  (4.12)  has  been  implemented 
for  the  data-adaptive  iterative  algorithm  to  address  this  problem.  This  modification 
involves  the  application  of  a  weighting  factor  to  the  current  update  of  Rn.  This 
weighting  factor  depends  on  the  index  of  iteration  n  such  that  its  value  approaches 
1  as  n  — ►  oo.  This  biasing  scheme  causes  the  influence  of  the  present  data  on  the 
parameter  estimates  to  increase  with  the  number  of  iterations.  The  modified  recursion 
is  given  by 

R^  =  AR»-i  +  (!L^-)  xnxj  .  (4.14) 

The  weighting  factor  on  the  second  term  in  (4.14)  is  added  to  provide  a  more  gradual 
shift  of  influence  to  recent  updates  than  would  be  possible  with  only  the  A  term 
present.  This  factor  is  not  as  important  when  a  highly  correlated  signal  is  processed. 

C.     SUMMARY 

The  adaptive  Toeplitz  approximation  algorithm  can  be  used  to  estimate  param- 
eters for  the  same  regions  of  support  as  in  the  fixed  data  case.  Additionally,  various 
schemes  for  pre- windowing  or  post-windowing  the  data  may  be  employed  to  improve 
the  estimates.  Other  factors  affecting  the  performance  of  the  algorithm  include  the 
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directional  scheme  for  moving  the  filter  over  the  data  array  and  the  criteria  used  to 
determine  convergence  of  the  parameter  estimates.  This  algorithm  has  been  used  for 
spectral  estimation,  image  noise  cancellation,  and  line  enhancement.  Experimental 
results  involving  these  applications  are  presented  in  the  next  two  chapters. 
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V.  2-D  AR  SPECTRAL  ESTIMATION 

A.  OVERVIEW 

One  application  where  the  use  of  two-dimensional  autoregressive  modeling  has 
proven  to  be  particularly  advantageous  is  that  of  spectral  estimation.  The  2-D  AR 
model  was  derived  in  Chapter  II.  This  model  is  used  to  find  the  2-D  AR  power 
spectral  estimate  Px(u)\,u>i)  which  is  given  by  [Ref.  4] 

A(«i,«i)  =  |J5T(fc*,u*)|aP„  ,  (5.1) 

where  Pw  is  the  power  spectral  density  of  the  input  and  H(u>i,uj2)  is  the  transfer 
function  of  the  2-D  AR  model.  The  input  is  white  noise  with  a  constant  power 
spectrum  of  amplitude  cr^.  Therefore,  we  can  write  (5.1)  as 

P*{UJU"2)  =  H  +  EE(,llfc2€.)a(^ib2)e-i(-^^)p  •  (5'2) 

The  key  to  finding  the  2-D  AR  power  spectral  estimate  is  determining  the  the  pa- 
rameters a(ki,k2).  The  next  section  provides  results  of  experimentation  with  the 
data-adaptive  iterative  algorithm  as  used  to  find  the  parameter  estimates  of  a  2-D 
signal  consisting  of  a  sum  of  sinusoids  in  additive  white  noise.  Results  using  the  di- 
rect inversion  method  and  fixed  data  iterative  Toeplitz  approximation  are  provided 
for  comparison. 

B.  EXPERIMENTAL  RESULTS 

The  performance  of  the  algorithms  discussed  in  the  previous  chapters  is  com- 
pared by  examining  the  estimated  spectral  densities  computed  from  the  parameters 
produced  by  each  method.  Each  algorithm  has  been  used  to  estimate  the  parameters 
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of  the  sinusoidal  input  array  generated  by 

x(nun2)    =    Ai  cos(27r/11n1  +  27r/12n2) 

+ A2  cos{2tt f2in1  +  2Trf22n2)  -\-  w{nun2)  (5.3) 

where  amplitudes  Ax  =  A2  =  y/2  and  /,_,  are  normalized  frequencies  in  the  range 
(0  <  fa  <  0.5).  The  frequencies  chosen  are  fn  =  /12  =  0.167  and  /2i  =  /22  =  0.333. 
Additionally,  the  data-adaptive  algorithm  is  tested  with  off-set  frequencies  to  evaluate 
its  peformance  for  that  case.  The  noise  term  iu(ni,  n2)  is  zero  mean  and  gaussian  with 
the  variance  <r£,  scaled  give  a  desired  signal-to-noise  ratio  (SNR).  The  SNR  (in  dB) 
is  defined  as  [Ref.  4] 

SNR=101og10£^-,  (5.4) 

,=1  Z(Tw 

where  N  is  the  number  of  sinusoids  present  and  A  is  the  amplitude  term.  The  SNR 
is  varied  by  holding  A,  constant  and  adjusting  the  value  of  a^.  The  algorithm's 
performance  using  a  3  x  3  filter  mask  is  evaluated  for  various  sizes  of  input  arrays 
and  the  common  regions  of  support,  including  combined-quadrant  (CQ).  The  results 
obtained  using  the  data-adaptive  algorithm  are  provided  for  the  pre-windowed  and 
covariance  methods.  The  data-adaptive  algorithm  is  tested  with  SNR's  of  10  dB  and 
0  dB.  It  should  be  noted  that  the  surface  plots  in  all  figures  have  been  rotated  90 
degrees  to  show  the  separation  of  the  spectral  peaks  more  clearly.  Contour  plots  are 
provided  to  better  judge  the  accuracy  of  the  estimates.  The  axes  of  the  contour  plots 
range  from  0-30  which  is  taken  from  a  32  x  32  frequency  domain  array  representing  a 
quadrant  of  a  64  x  64  point  2-D  FFT  output.  The  32  x  32  array  covers  the  frequency 
range  from  0  to  7r. 

1.      Estimates  by  Direct  Inversion 

The  first  case  examined  is  that  involving  direct  inversion  of  the  autocorre- 
lation matrix  to  find  the  least  squares  AR  parameters  from  the  normal  equations.  An 
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8x8  input  array  was  used  with  l*1  quadrant  quarter-plane  support.   The  resulting 
spectral  density  estimate  is  shown  in  Figure  5.1.   This  plot  shows  spectral  peaks  at 


E      15 


Fl 
Figure  5.1:  Spectral  estimate  using  direct  matrix  inversion;  SNR=10  dB. 

21.3  and  10.7  on  the  Fl  and  F2  axis  respectively,  which  correspond  to  the  normalized 
frquencies  fn  =  /12  =  0.167  and  f2\  =  J27  =  0.333  of  the  input  data.  The  true 
locations  of  the  frequencies  are  indicated  by  crosses. 

2.      Spectral  Estimates  using  the  Iterative  Method  for  Fixed  Data 
a.      Quarter-Plane  Support 

Figure  5.2  shows  the  spectral  density  of  a.-(ni,n2)  as  estimated  using 
the  fixed  data  iterative  method  with  first  quadrant  QP  support.  As  in  the  previous 
case  the  frequencies  estimated  are  very  close  to  the  true  frequencies.  It  is  interesting 
to  note  however,  that  contrary  to  what  might  be  expected,  the  quality  of  the  estimate 
appears  to  degrade  as  the  number  of  iterations  is  increased.  The  estimate  produced 
after  one  iteration  is  clearly  superior  to  the  estimate  obtained  after  ten  iterations. 
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Tliis  phenomenon  is  consistent  with  that  experienced  in  the  1-D  case  where  the 
iterative  solution  would  begin  to  diverge  after  a  certain  number  of  iterations  in  cases 
with  small  data  sets  [Ref.  12]. 
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Figure  5.2:  Fixed  Data,  1st  Quadrant  QP  Support,  SNR=10  dB. 

Results  using  second  quadrant  QP  support  proved  to  be  generally 
better  for  this  placement  of  the  sinusoids  than  those  obtained  from  the  first  quadrant. 
An  example  is  provided  in  Figure  5.3. 

b.      Combined  Quadrant  Support 

Spectral  density  estimates  using  QP  support  demonstrate  a  tendency 
to  distort  elliptically.  This  distortion  is  evident  in  the  first  quadrant  QP  case  shown 
in  Figure  5.2.    A  method  to  compensate  for  this  elongation  of  the  spectral  peaks 
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Figure  5.3:  Fixed  Data,  2nd  Quadrant  QP  Support,  SNR=10  dB. 
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involves  combining  the  first  and  second  quadrant  estimates  [Ref.  13]  to  give 

2a.2. 


P,= 


(5.5) 


The  plots  in  Figure  5.4  indicate  that  this  method  yields  accurate  estimates.  However, 
despite  some  improvement,  the  elliptical  elongation  at  the  base  of  the  spectral  peaks 
resulting  from  the  influence  of  the  first  quadrant  estimate  is  still  quite  noticeable. 
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Figure  5.4:  Fixed  Data,  Combined  Quadrant  Support,  SNR=10  dB. 

3.      Spectral  Estimates  using  the  Data-Adaptive  Iterative  Method 
a.      Pre- Windowed  Input  Arrays 

The  input  arrays  were  pre-windowed  by  adding  Px  —  1  rows  and 
columns  of  zeros  to  the  top  and  left  side  of  the  array,  where  P\  is  the  dimension 
of  the  filter  mask.  The  purpose  of  pre-windowing  the  data  is  to  compute  the  pa- 
rameter estimates  using  all  the  input  data.   Figure  5.5  shows  the  spectral  estimates 
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for  1st,  2nd  and  CQ  support  with  an  8  x  8  input  array.  In  Figure  5.6  we  present  the 
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Figure  5.5:  8x8  input  array  (pre-windowed),  SNR=10  dB. 

spectral  estimates  when  only  16  data  points  (4x4  array)  were  available.  The  spectral 
estimates  obtained  from  the  8x8  input  array  were  quite  acceptable,  but  some  bias 
is  evident  in  the  first  quadrant  estimate  which  is  carried  over  to  the  CQ  estimate. 
The  estimates  obtained  from  the  4x4  array  demonstrate  the  value  of  computing  the 
CQ  estimate.  In  this  case  the  true  signal  frequencies  were  unrecognizable  until  the 
combined  quadrant  estimate  was  obtained. 

A  final  test  of  the  pre-windowed  algorithm  consisted  of  running  the 
filter  mask  only  half  way  through  an  8  x  8  input.  The  estimates  are  computed  from 
the  upper  4x8  section  of  the  observed  data.  The  result  of  this  experiment  is  given  in 
Figure  5.7.  These  results  are  similar  to  that  of  the  4x4  case  and  serve  to  reinforce  the 
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Figure  5.6:  4x4  input  array  (pre-windowed),  SNR=10  dB. 
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Figure  5.7:  4x8  portion  of  the  input  array  (pre-windowed),  SNR=10  dB. 
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importance  of  the  Combined  Quadrant  technique  for  2-D  AR  spectrum  estimation. 
b.      Estimates  using  the  Covariance  Method 

As  explained  in  section  B  of  Chapter  II  the  covariance  method  uses 
only  the  observed  data.  As  a  result,  Pi  —  1  rows  and  columns  of  the  input  array  are 
not  available  when  computing  the  parameter  estimates.  However,  the  experimental 
results  obtained  in  this  work  indicate  that  the  covariance  method  achieves  better 
spectral  peak  resolution.  Additionally,  less  bias  is  observed  in  the  spectral  estimates. 
This  characteristic  of  the  covariance  method  is  evident  in  Fig  5.8.    In  this  case  the 
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Figure  5.8:  8x8  input  array  (Covariance  Method),  SNR=10  dB. 

best  spectral  estimate  is  clearly  provided  by  the  Combined  Quadrant  method.  The 
bias  present  in  the  pre-windowed  first  quadrant  support  example  has  been  completely 
eliminated  by  using  the  covariance  method. 
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The  next  example  demonstrates  a  unique  characteristic  of  the  data- 
adaptive  iterative  Toeplitz  approximation  algorithm.  In  this  example  the  spectral 
estimate  is  computed  using  a  3  x  3  filter  mask  with  a  4  x  4  input  array.  This  means 
that  only  four  elements  of  the  input  data  are  available  to  the  filter  to  determine 
the  parameter  estimates.  This  result  is  of  particular  interest  because  the  correlation 
matrix  is  actually  rank  deficient,  and  the  problem  cannot  be  solved  by  direct  inversion. 
Nevertheless,  the  iterative  Toeplitz  approximation  method  does  provide  a  means  to 
compute  the  estimate.    The  results  are  given  in  Figure  5.9.    Although  the  spectral 
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Figure  5.9:  4x4  input  array  (Covariance  Method),  SNR=10  dB. 

estimates  for  the  first  quadrant  and  Combined  Quadrants  are  somewhat  weak,  the 
second  quadrant  estimate  exhibits  remarkable  resolution  and  accuracy. 

As  a  further  test  of  the  algorithm,  a  signal  with  ofFset  frequencies 
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/n  =  0.250,  /12  =  0.125,  f2\  =  0.300,  and  /22  =  0.400  is  generated  and  processed. 
The  spectral  estimates  resulting  from  the  computed  model  parameters  are  plotted 
in  Figure  5.10.     For  this  case,  the  best  spectral  estimate  was  obtained  using  first 
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Figure  5.10:    8x8  input  array,  offset  frequencies,  (Covariance  Method), 
SNR=10  dB. 

quadrant  suppport. 

The  spectral  estimates  for  a  signal-to-noise  ratio  of  0  dB  is  ploted 
in  Figure  5.11.  The  accuracy  of  the  estimates  is  noticeably  degraded  in  this  case. 
However,  a  reasonable  indication  of  the  signaj  frequency  content  can  be  seen  when 
the  combined-quadrant  method  is  used. 

The  final  case  using  NSIIP  support  is  included  for  completeness.  The 
filter  mask  had  dimensions  P\  =  P2  =  4  and  L2  =  —2.  The  parameters  were  estimated 
for  a  16  x  16  input  array.  The  NSHP  spectral  estimate  is  plotted  in  Figure  5.12.  The 
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Figure  5.11:  8x8  input  array  (Covariance  Method),  SNR=0dB. 
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Figure  5.12:  16  x  16  input  array  (Covariance  Method),  SNR=10  dB;  non- 
symmetric  half-plane  support. 


40 


resulting  spectral  estimate  shown  in  Figure  5.12  is  quite  accurate,  but  this  comes 
at  the  cost  of  greater  complexity  for  the  filter.  NSHP  support  may  be  used  as  an 
alternative  to  combined- quadrant  support  in  some  cases. 

C.      SUMMARY 

The  data-adaptive  iterative  Toeplitz  approximation  algorithm  has  proven  to 
be  a  viable  means  for  2-D  AR  spectrum  estimation.  Numerous  factors  must  be 
considered  when  using  this  algorithm.  In  particular,  the  filter  mask  size  should  be 
chosen  to  correspond  to  the  model  order  of  the  signal  being  processed.  Although 
our  overall  experience  has  found  that  that  CQ  support  produced  the  best  results, 
no  region  of  support  seemed  to  have  a  clear  advantage  over  the  others  for  all  cases. 
This  indicates  that  the  optimum  region  of  support  must  be  evaluated  for  the  specific 
application  that  the  algorithm  is  implemented  for. 
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VI.  RESTORATION  OF  IMAGES 

A.      OVERVIEW 

The  characteristics  of  typical  images  vary  greatly  from  one  region  of  the  image 
to  the  next.  For  example,  a  region  of  an  image  containing  a  crowd  or  a  building  may 
have  detailed  variations  in  intensity,  while  another  region  representing  the  sky  in  the 
background  will  be  essentially  uniform  in  intensity.  Two-dimensional  data-adaptive 
filtering  is  particularly  suited  for  processing  data  of  this  nature.  The  idea  of  adapting 
the  processing  to  the  local  characteristics  of  an  image  is  advantageous  for  many  image 
processing  applications,  including  image  enhancement  and  restoration. 

There  are  two  approaches  to  adaptive  signal  processing.  The  first  approach 
involves  adapting  the  iilter  output  for  each  pixel.  This  is  known  as  pixel-by-pixel 
processing  [Ref.  4].  In  this  scheme,  the  processing  method  is  based  on  the  local 
characteristics  of  the  image,  degradation,  and  any  other  pertinent  information  con- 
tained in  the  pixel's  neighborhood  region.  This  approach  offers  the  greatest  degree 
of  fiexiblity  to  the  adaptive  process,  but  has  the  highest  computational  complexity. 

When  a  more  computationally  efficient  method  is  desired,  a  second  approach 
known  as  subimage-by-subimage  or  block-by-block  processing  can  be  used  [Ref.  4].  In 
this  approach,  the  image  array  is  divided  into  subimages  and  space-invariant  tech- 
niques are  used  to  process  the  data.  This  method  can  result  in  some  discontinuities 
being  present  in  the  processed  image.  The  size  of  the  subimages  directly  affects  the 
quality  of  the  output  image.  The  block-by-block  method  provides  a  faster  means  of 
processing,  however,  and  may  be  desirable  in  some  applications. 

In  the  next  two  sections  of  this  chapter  the  data-adaptive  Toeplitz  approxima- 
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tion  algorithm  is  implemented  in  image  noise  cancellation  and  line  enhancer  modes. 
For  both  of  these  applications  the  algorithm  has  been  modified  to  perform  pixel-by- 
pixel  processing. 

B.     IMAGE  NOISE  CANCELLATION 

The  usual  method  of  finding  the  estimate  of  a  signal  in  noise  is  to  pass  the 
corrupted  signal  through  a  system  that  serves  to  suppress  the  noise  while  leaving 
the  desired  signal  relatively  unchanged.  The  noise  canceler  developed  by  Widrow 
[Ref.  14]  is  an  example  of  such  a  system.  Adaptive  noise  cancellation  is  a  variation 
of  optimal  filtering  that  can  be  used  for  image  restoration  to  some  advantage.  In 
particular,  when  a  recieved  image  is  obscured  by  noise  or  another  image  transmitted 
from  another  location  as  shown  in  Figure  6.1,  an  adaptive  noise  canceler  can  be 
employed  to  extract  the  desired  image  from  the  composite  signal.  The  noise  canceler 

unwanted   transmission 


reference   sensor 


desired   transmission 


PROCESSING 
UNIT 
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Figure  6.1:   A  possible  scenario  for  application  of  a  noise  canceler  as  used 
for  image  restoration. 


makes  use  of  an  auxiliary  or  reference  input  derived  from  one  or  more  sensors  located 
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at  points  in  the  noise  field  where  the  desired  signal  is  weak  or  undetectable.  This 
input  is  filtered  and  subtracted  from  the  primary  input  containing  both  desired  and 
unwanted  signals.  A  block  diagram  of  the  noise  canceler  is  shown  in  Figure  6.2.  The 
reference  signal  and  the  undesired  part  of  the  primary  input  are  correlated.  For  this 
reason,  the  unwanted  signal  or  noise  is  eliminated  or  attenuated  by  cancellation. 
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Figure  6.2:  Adaptive  Noise  Canceler  Block  Diagram. 

To  evaluate  the  performance  of  the  data-adaptive  iterative  Toeplitz  approxi- 
mation algorithm  used  in  a  noise  canceler  configuration,  a  256  x  256  composite  im- 
age array  was  created  from  two  distinct  images.  This  corrupted  image  is  shown  in 
Figure  6.3.  The  algorithm  described  in  Chapter  IV  was  modified  to  include  a  cross- 
correlation  term  given  by 

rn  =  Arn_i  +  I J  dnxn  (6.1) 

where  dn  is  the  element  of  the  composite  signal  corresponding  to  the  nth  element  of 
the  reference  signal  being  processed  by  the  filter.  This  cross-correlation  term  is  used 
to  form  the  initial  parameter  estimate  for  each  pixel  being  processed. 
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Figure  6.3:    Corrupted  image  array  provided  as  primary  input  to  noise 
canceler. 

The  resulting  error  e(ni,ri2)  computed  at  the  output  of  the  system  becomes  the 
restored  image.  This  relationship  is  given  by 

e{nun2)  =  [d{nun2)  +  u{nun2)]  -  u(ni,n2)  ,  (6.2) 

where  d  is  the  desired  image,  u  is  the  unwanted  image  or  noise  and  u  is  the  filtered 
estimate  of  u. 

In  the  simulation  the  noise  signal  added  to  the  desired  image  is  scaled  and 
lowpass  filtered  to  simulate  the  effect  of  having  the  primary  and  reference  sensors 
physically  separated.  The  output  of  the  noise  canceler  is  shown  in  Figure  6.4.  The 
unwanted  image  has  been  noticeably  attenuated  and  the  desired  image  is  clearly 
visible.  The  original  images  are  provided  in  Figure  6.5  as  a  reference  for  evaluation 
of  the  algorithm's  performance. 

C.      ADAPTIVE  LINE  ENHANCER 

A  special  case  of  adaptive  noise  canceling  occurs  when  only  the  signal  that  is 
corrupted  by  noise  is  available.  In  this  case  the  recieved  signal  is  delayed  and  used  as 
the  reference  signal.  The  reference  signal  may  be  expressed  as  y(nu  n2)  =  x(rii—  6,  n2), 
where  8  is  the  amount  of  spatial  delay  used.    The  block  diagram  for  this  system  is 
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Figure  6.4:  Processed  image  array  after  noise  cancellation. 


Figure  6.5:  The  original  images  -  undesired  and  desired. 
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shown  in  Figure  6.6.   The  main  function  of  the  delay  parameter  6  is  to  remove  any 
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Figure  6.6:  The  adaptive  line  enhancer  block  diagram. 

correlation  that  may  exist  between  the  noise  component  in  the  original  input  signal 
and  the  noise  component  of  the  delayed  adaptive  filter  input  [Ref.  10].  The  filter  will 
respond  by  cancelling  any  components  of  the  main  signal  x(ni,n2)  that  are  in  any 
way  correlated  with  the  secondary  signal  y{n\,n2)  [Ref.  11].  The  remaining  noise 
component  at  the  output  of  the  filter  is  then  subtracted  from  main  signal  resulting 
in  the  removal  of  the  additive  noise  in  the  signal.  In  general,  the  delay  6  should 
be  chosen  such  that  it  is  approximately  half  the  length  of  the  correlation  sequence 
[Ref.  11].  The  input  image  and  output  of  the  adaptive  line  enhancer  are  shown  in 
Figure  6.7. 

D.      SUMMARY 

In  this  chapter,  it  has  been  demonstrated  that  the  data  adaptive  Toeplitz  ap- 
proximation algorithm  may  be  used  in  image  processing  applications.  While  pixel-by- 
pixel  processing  was  used  in  the  above  examples,  it  may  be  desirable  to  combine  this 
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Figure  6.7:  Input  and  output  of  the  data-adaptive  line  enhancer. 

approach  with  block-by-block  processing  to  obtain  better  results.  This  would  entail 
dividing  the  image  into  smaller  blocks  before  processing  and  then  recombining  the 
processed  blocks.  Additionally,  reprocessing  of  the  array  estimates  may  be  a  viable 
means  of  improving  image  quality. 
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VII.  CONCLUSIONS 

The  data-adaptive  iterative  Toeplitz  approximation  algorithm  is  readily  im- 
plemented from  the  iterative  methods  previously  used  for  fixed  data-arrays.  The 
experimental  results  indicate  that  this  algorithm  is  well-suited  for  a  wide  variety  of 
2-D  signal  processing  applications.  Practically  speaking,  the  iterative  method  can 
be  successfully  implemented  for  any  application  where  the  use  of  Wiener  filter  based 
adaptive  filters  has  been  successful.  While  more  extensive  evaluation  is  necessary 
before  this  method  can  be  said  to  have  a  clear  advantage  over  other  methods,  there 
is  enough  evidence  of  its  capabilities  to  warrant  further  investigation  in  this  area. 

A.     PERFORMANCE  EVALUATION  SUMMARY 

As  stated  earlier,  autoregressive  model  parameters  may  used  to  estimate  the 
spectral  content  of  2-D  arrays  with  high  resolution.  In  experiments  the  performance 
of  the  data-adaptive  iterative  algorithm  matched  that  of  the  direct  inversion  method, 
and  in  some  cases,  improved  upon  the  performance  of  the  fixed  data  iterative  method. 
This  improvement  may  be  due  to  the  way  in  which  the  iteration  is  carried  out  in  the 
adaptive  case  versus  the  fixed  case.  In  the  fixed  case  the  iteration  is  carried  out  using 
a  correlation  matrix  formed  from  all  the  data  at  once.  In  the  data-adaptive  case  the 
correlation  matrix  is  formed  recursively  with  new  data  incorporated  at  each  itera- 
tion. Of  particular  note,  was  the  ability  of  the  data-adaptive  algorithm  to  estimate 
frequencies  using  very  small  data  sets  and  its  ability  to  produce  an  estimate  even 
when  the  correlation  matrix  is  not  of  full  rank  .  Additionally,  its  performance  in  a 
high  noise  environment  was  noteworthy. 

The  algorithm  also  proved  to  be  a  viable  method  for  data-adaptive  image 
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restoration.  The  nature  of  this  application  differs  greatly  from  that  of  spectral  estima- 
tion in  that  the  data  sets  are  much  larger  and  generally  less  correlated.  Encouraging 
results  were  obtained  for  the  data-adaptive  noise  canceler  and  line  enhancer  with  only 
minor  modifications  of  the  same  algorithm  that  was  used  for  the  spectrum  estimates. 

B.     RECOMMENDATIONS  FOR  FUTURE  STUDY 

The  primary  objective  of  this  thesis  was  to  develop  the  data-adaptive  iterative 
Toeplitz  algorithm  and  obtain  experimental  results  to  evaluate  its  potential  for  use  in 
applications  that  involve  2-D  AR  parameter  estimation.  There  are  many  aspects  re- 
garding this  work  that  require  more  in-depth  study  and  several  applications  that  may 
benefit  in  the  use  of  this  algorithm.  In  particular,  more  detailed  analysis  of  model 
orders  and  their  relation  to  the  algorithm's  ability  to  form  spectral  estimates  from 
very  small  input  arrays  needs  to  be  carried  out.  The  use  of  singular  value  decompo- 
sition (SVD)  and  eigenvalue  analysis  of  the  observed  data  may  provide  some  insights 
regarding  the  performance  of  the  iterative  algorithm.  Also,  further  investigation  as  to 
why  one  region  of  support  (ROS)  provides  greater  resolution  than  another,  may  yield 
information  that  would  permit  the  choice  of  an  optimum  ROS  prior  to  processing 
the  data.  Further  analysis  for  choosing  values  of  the  forgetting  factor  A  and  how  it 
relates  to  the  statistical  nature  of  the  data,  is  desirable  as  well.  This  may  lead  to  a 
scheme  for  on-line  adjustment  of  A,  which  would  be  particularly  advantageous  when 
processing  images.  Other  future  work  could  explore  the  use  of  different  data  win- 
dowing schemes  and  data-ordering  schemes  for  passing  the  filter  mask  over  the  data. 
One  possible  method  that  has  not  been  investigated  is  an  expanding  square  starting 
at  one  corner  of  the  array.  This  method  seems  well  suited  for  block-by-block  image 
processing.  Finally,  the  data-adaptive  iterative  Toeplitz  approximation  algorithm  has 
the  potential  for  use  in  non-linear  applications,  such  as  the  identification  of  Volterra 
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